function trans_correct(dfile, n)
%remove the transient from traces
% global CONTROL update 3/23/02 PBM

global CURRENT VOLTAGE WCHAN
global CONTROL
%clean_current=current_in; % first, just make output same as input
QueMessage(sprintf('Transient Corrections'));
[recs,pts]=size(CURRENT);

[trans_templ]=find_transient(dfile, n);
[v_steps]=find_tr_steps;
if(isempty(v_steps) | isempty(trans_templ)) % did'nt find any?
	QueMessage(sprintf('Transient Corrections - no steps or template'));
	return
end;
	
%clean_up
%now convolve transient template with voltage steps to get current
for i=1:recs
   transient(i,:)=conv(v_steps(i,:),trans_templ);
   transient(i,pts)=0;
   CURRENT(i,1:pts)=CURRENT(i,1:pts)-transient(i,1:pts); %shift by one
end
return

% *****************************************************************
% 
% *****************************************************************
function [v_steps]=find_tr_steps()
%finds the voltage steps in the command and outputs a dirac sequence

global VOLTAGE

v_steps=[]; % default return value is NULL..
v_steps1=diff(VOLTAGE')';
[v_stpeaksr,v_stpeaksc]=find(abs(v_steps1) > 1);
[vr,vc]=size(v_steps1);
v_steps=zeros(vr,vc);
if(length(v_stpeaksr) > 200 | length(v_stpeaksc) > 200)
   return;
end;

v_steps(v_stpeaksr,v_stpeaksc)=v_steps1(v_stpeaksr,v_stpeaksc);

return

% *****************************************************************
% 
% *****************************************************************
function [trans_templ]=find_transient(dfile, n)
%find and make the transient template
%length of transient template: 6 ms => fill up to nearest
%32/64 etc with zeros
% modified from pok version, by pbm, 10/13/98

global CURRENT VOLTAGE
global CONTROL

trans_templ=[]; % default return value is NULL
step = convert_steps_to_times(n);
rate = dfile.rate(1)*dfile.nr_channel/1000;
tr_lengthms=number_arg(CONTROL(n).transwdur);
tr_length=floor(tr_lengthms./rate);

trwin = number_arg(CONTROL(n).transwin);
tr_start=step(trwin).begin;
tr_end = tr_start+number_arg(CONTROL(n).transwdur);
tr_length = floor((tr_end - tr_start)./rate);
QueMessage(sprintf('Transient Correction in window: %8.2f - %8.2f, tr_length = %5d', tr_start, tr_end, tr_length));
if tr_end-tr_start < tr_lengthms
   tr_end=tr_start+tr_lengthms+5;
end

%find steps in [tr_start,tr_end], when v is in [v1, v2]
t_start=floor(tr_start./rate);
if(t_start<=0)
t_start=1;
end
t_end=floor(tr_end./rate);
if(t_end > size(VOLTAGE, 2)-1)
    t_end = size(VOLTAGE, 2)-1;
end;
%disp(sprintf('rate: %f start: %d  end: %d', rate, t_start, t_end ))
%average voltage in tr_start tr_end

lwin = number_arg(CONTROL(n).leakwin);
lt1 = step(lwin).begin+0.5*step(lwin).dur;
lt2 = step(lwin).begin+step(lwin).dur*0.95;
lt3 = 0.5*step(1).dur; % use half of the holding step duration
t1=floor(lt1./rate); %begin of mean
t2=floor(lt2./rate); %end of mean
if(t1<=0) t1=1; end;
if(t2<=0) t2=1; end;

v1=number_arg(CONTROL(n).vlmax); %usually -90 mV
v2=number_arg(CONTROL(n).vlmin); %usually -65 mV
if(v1 > v2)  % swap arguments
   t = v2;
   v2=v1;
   v1=t;
end
[recs pts] = size(VOLTAGE);
if(t1 > pts | t2 > pts)
   QueMessage(' Bad Window for Leak Correction');
   return;
end;
volt=mean(VOLTAGE(:,t1:t2)');
%figure out which traces in the voltage window
v_ind=find(volt > v1 & volt < v2);
if length(v_ind) < 2
   %no data in range
   p(1)=0;
   p(2)=0;
   tra=0;	
   QueMessage('No data to correct!');
   return;
end
av_volt=mean(VOLTAGE(v_ind,t_start:t_end+1));
av_curr=mean(CURRENT(v_ind,t_start:t_end+1));
t_length=t_end-t_start+1;
%find bginning of step
dvolt=diff(av_volt);
[mdv,mdi]=min(dvolt);  %step occured at mdi
% plot (av_volt);
%now get current from mdi to tr_end
if mdi+tr_length > t_end
   tr_length=t_end-mdi;
end

clear trans_templ
trans_templ(1)=0;  %set beginning to zero
%disp(sprintf('mdi: %d, tr_length: %d, ', mdi, tr_length))
in_size=length(av_curr)
if(tr_length-1 > in_size)
   trans_templ=zeros(1, max(tr_length));
   QueMessage('Check windows for transient correction');
   QueMessage('Correction array set to zeros (no correction applied)');
else
   trans_templ(1:tr_length)=av_curr(1:tr_length);
end;

%scale by the voltage step
trans_templ=trans_templ./mdv;

%now pad with zeros
ttl=(floor(length(trans_templ)/16)+1)*16;
trans_templ(tr_length+2:ttl)=0;

return
